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Abstract 

An improved 4-node quadrilateral assumed-stress hybrid shell element with drilling 
degrees of freedom is presented. The formulation is based on Hellinger-Reissner varia- 
tional principle and the shape functions are formulated directly for the 4- node element. 
The element has 12 membrane degrees of freedom and 12 bending degrees of freedom. 
It has 9 independent stress parameters to describe the membrane stress resultant field 
and 13 independent stress parameters to describe the moment and transverse shear 
stress resultant field. The formulation encompasses linear stress, linear buckling, and 
linear free vibration problems. The element is validated with standard test cases and 
is shown to be robust. Numerical results are presented for linear stress, buckling, and 
free vibration analyses. 
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Introduction 


Today, the literature abounds with new shell element formulations and the key is- 
sue is the robustness of individual elements. Factors that influence the robustness 
of elements are spurious zero-energy modes, locking, convergence characteristics, ele- 
ment invariance, and sensitivity to mesh distortion. Considerable progress has been 
made over the years in addressing some of these factors and in developing suitable 
assessment procedures and test cases to evaluate and validate the new elements. In 
the process, several new techniques and new element formulations such as assumed- 
stress, reduced-integration, incompatible modes, and assumed-strain have evolved, 
that successfully address some of the problems associated with shell elements. 

On a geometrical aspect, the shell element formulations that are currently avail- 
able in literature can be broadly classified into three categories: curved elements based 
on classical shell theories, degenerated solid shell elements, and flat shell elements. 
The faceted representation of the shell geometry using flat shell elements is perhaps 
the simplest of these approaches. These elements combine the membrane and bend- 
ing behavior of plate elements. While the approximate geometrical representation is 
immediately evident, particularly in 4-node bilinear elements, the simplicity of the 
formulation, the convergence characteristics, and full rank has made this approach 
very attractive. The inclusion of transverse shear effects with the aid of Reissner- 
Mindlin kinematics, and more recently, the inclusion of drilling degrees of freedom 
[1, 2], have significantly improved the performance of flat shell elements. 

The assumed-stress hybrid formulations pioneered by Pian and his coworkers 
[3] is variationally consistent and has led to the development of several powerful 
element formulations. In general, the mixed/hybrid formulations are computationally 
more intensive than the displacement-based formulations; however, using symbolic 
manipulation (e.g., see [4]), and exploiting the banded structure of matrices, the 
computational cost could be reduced significantly. 
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The objective of this paper is to present the formulation of an improved 4-node 
assumed-stress hybrid shell element with drilling degrees of freedom. The formula- 
tion encompasses linear stress, linear buckling, and linear free vibration problems. 
The element is validated using standard test cases as well as a few other interesting 
problems for linear stress, buckling, and free vibration analyses. Wherever possible, 
relevant matrices are generated symbolically, and this has helped save considerable 
computational time. However, the details of the symbolic manipulation will be dealt 
with in a subsequent paper. 

Drilling Degrees of Freedom 

Drilling degrees of freedom are defined as the rotational degrees of freedom normal 
to the plane of the element. The drilling rotation (0 Z ) is not represented in the mem- 
brane kinematics of the shell. Consequently, the drilling degree of freedom is not 
included at the element level; however at the global or assembled level, the degrees 
of freedom associated with the structure may include the drilling degrees of freedom. 
While assembling the element stiffness matrices to form the global stiffness matrix, if 
the elements connected to a particular node happen to be coplanar, then the drilling 
rotation at that node is not resisted. This leads to a singularity in the stiffness matrix 

[5] . This singularity could clearly be avoided if the drilling degree of freedom is in- 
cluded in the variational statement or in the finite element approximations as will be 
discussed subsequently. Yet another motivation to include the drilling degree of free- 
dom stems from difficulties in modeling stiffened panels, folded plate structures, etc. 
The lack of drilling rotation in such instances results in an inadequate representation 
of the structural response. 

Generally, the drilling degrees of freedom are suppressed at the beginning of the 
analysis since they do not enter the kinematic description of the problem. Zienkiewicz 

[6] associates a fictitious couple M z with the rotation at the element level to resolve 
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the singularity problem. While this does take care of the singularity in the stiff- 
ness matrix, it does not utilize the drilling degrees of freedom to improve the finite 
element approximations. However, there have been numerous attempts in formulat- 
ing elements with drilling degrees of freedom, and two significant approaches have 
emerged successful. One approach is to introduce the drilling rotation in the varia- 
tional statement as an independent variable, while the other approach is to introduce 
the drilling rotation in the approximations of in-plane displacements. 

Drilling Rotation in the Variational Statement 

Reissner [7] presented a mixed variational principle in which the definition of true 
rotations as the skew-symmetric part of the displacement gradients are relaxed and 
later enforced in the variational statement as constraints through the introduction of 
Lagrange multipliers. On imposing the stationary conditions on the functional, the 
Lagrange multipliers are identified as the skew-symmetric part of the stress tensor. 
This formulation introduced a rotation field independent of displacements. Hughes 
and Brezzi [8] modified Reissner’s functional by including an additional term that 
stabilizes the functioned in discrete approximations. Ibrahimbegovic et al. [9, 10] 
used the modified formulation of Hughes and Brezzi and developed quadrilateral 
membrane and flat shell elements. They used independent interpolation fields for 
rotations and Allman-type [1] shape functions for the in-plane displacements. Iura 
and Atluri [11] have developed an element where the rotations are interpolated from 
true nodal rotations evaluated from the skew-symmetric part of the displacement 
gradient at the nodes. 

Drilling Rotations in Finite Element Approximations 

In early attempts, the drilling rotation was introduced in the shape functions 
using a cubic displacement function (e.g., see Robinson [12]). However, Irons and 
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Ahmad [13] pointed out that this method had some difficulties and convergence was 
not assured. The elements formed in such a manner force the shearing strain to be 
zero at the nodes and do not pass the patch test. Recently, the drilling rotations have 
been introduced in the finite element approximations successfully using quadratic dis- 
placement functions [1,2,14-18]. In this approach, the membrane element is internally 
assumed to be an 8-node element (four corner nodes and four midside nodes) with 16 
degrees of freedom. The stiffness matrix of this 8-node element is then condensed to 
that of a 4-node element (four corner nodes) with 12 degrees of freedom by associat- 
ing the displacement degrees of freedom at the midside nodes with the displacement 
and rotational degrees of freedom at the corner nodes. In effect, the in-plane dis- 
placements are dependent on in-plane rotations. MacNeal and Harder [16] have used 
the above approach to develop a 4-node membrane element with selectively reduced 
integration. Yunus et al. [17] have also used this approach to develop assumed-stress 
hybrid elements. Aminpour [18, 19] has developed assumed-stress hybrid shell ele- 
ments using the Allman-type shape function construction. The formulation of the 
first element is based on Hellinger-Reissner variational principle. The formulation 
of the second element is based on modified complementary energy principle and the 
stress field is expanded in Cartesian coordinates. 

Assumed-Stress Hybrid Formulation 

The element developed in this paper is based on Hellinger-Reissner variational prin- 
ciple. The ensuing discussion pertains to a solid elastic continuum V, with boundary 
5. Let S„ be the section of the boundary where tractions are prescribed; let S u be the 
section of the boundary where displacements are prescribed. The Hellinger-Reissner 
functional can be written as 

n„ B = -[ / M r [D]M dv + f M t [£|M dv - [ {u) T {Q ds, (i) 

2 Jv JV Js a 
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where [ D ] is the compliance matrix, \C\ is the linear differential operator on the 

displacements {u} to produce strains, and {£„} is the prescribed tractions on S„. 

The approximations for stresses and displacements can now be incorporated in 
the functional. The stress field is described in the interior of the element as 

m = [pm, ( 2 ) 

and a compatible displacement field is described by 

M = [N}{qh (3) 

where [P] and [N] are matrices of stress and displacement interpolation functions, 
respectively, and {/?} and {< 7 } are the unknown stress and nodal displacement param- 
eters, respectively. Substituting the stress and displacement approximations [Equa- 
tions ( 2 ) and (3)] in the functional Iljfii [Equation ( 1 )] results in 

n™ = -\m T \»m + - { ? } t {f>, (4) 

where 

[H] = j v [P} T [D}[P\ dV , 

m = j v \p\ T ww w , 

{F} = / [JV] T {( 0 } dS . (5) 

JS* 

Now imposing stationary conditions on the functional with respect to the stress pa- 
rameters {/?} gives 

m = [ffr'mw- (6) 

On substituting Equation ( 6 ) into Equation (4), the functional reduces to 

n HB = i{9} r [JT]{ 9 } - M r {F}, (7) 

where [if], the elemental stiffness matrix, is given by 

[ifl = (TflflT'lT]. (8) 
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Imposing the stationary conditions on the functional with respect to nodal displace- 
ment unknowns {7}, results in a system of linear algebraic equations of the form 

[km = m- (9) 

Finite Element Approximations 

Displacement Field 

In the present paper, the Allman-type shape functions are used. However, the 
formulation is direct rather than forming the stiffness matrix of an 8-node element 
and then condensing it to a 4-node element stiffness matrix. All four edges of the 
quadrilateral may be treated as 2-D beam elements with shear deformation, and the 
general interpolation function for any typical edge can be expressed as in reference 
[18] by 

= 5(1 - (K + |(i + 0«, + - 0) (*., - *-). 

1 1 A 2 * * 

»° = 5(1-0* + 5(1 + 0-, - -^(1 -<’)(*«- *») (10) 

and 

Ax» = Xj - Xi , 

Ay; = yj - Vi , 

where i = 1,2, 3, 4 and j = 2,3,4, 1 in that order. As mentioned earlier, the drilling 
rotations are not interpolated independently. The out-of-plane displacement and out- 
of-plane rotations are given by 

= 5 O -O'* + 5(1 + 0-, - ^(1 

A T . 

+ -^(1 - 0) («„ - <*) 

». = j(l-0«- + 5(1+0*.,. 

*, = 5(1 - 0*- + 5(1 + 0 «.,• (11) 
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Note that till three displacement degrees of freedom and two in-plane rotations 
are approximated by quadratic and linear functions, respectively. This difference in 
the order of approximating polynomials between displacements and rotations auto- 
matically alleviates shear locking. Equations (10) and (11) can be extended directly 
to all four sides of the quadrilateral element to obtain the interpolation functions 
for the element. The generalized displacements represented in terms of interpolation 
functions are 




Ni = 

^(1+60(1 + ViV)i i = 1,2, 3, 4 

(13) 



| JU-f’JO + iCT) ; * = i,3 



n: = 

1 1(1 -*,’)(!+«) ; * = 2,4 

(14) 

and 


i + 1 ; i = 1, 2, 3 




j = < 

(15) 



1 ; i = 4 



In the above set of equations £ and 7 / are the element natural coordinates and (i and 


rji are the values of £ and 77 at node i. Note that the drilling degrees of freedom (0 X ;) 
are not true nodal rotations but are referred to as “rotational connectors” [ 1 ] and 
enter the finite element approximations as differences in rotational connectors. 

As opposed to conventional 4-node shell elements which have 20 degrees of free- 
dom, elements with drilling degrees of freedom have 24 degrees of freedom. The 
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membrane part of the element has 12 degrees of freedom. Associated with these 12 
degrees of freedom are 12 independent deformation modes. Three of these modes 
are rigid-body modes and three other modes represent constant strain states. Of the 
remaining six modes, five represent higher-order strain states and the last mode is 
a spurious zero-energy mode. The spurious zero-energy mode associated with this 
formulation is caused by the drilling rotations or rather the rotational connectors. 
Note from Equations (10), (11) and (12) that the differences in rotations, and not the 
drilling rotations themselves, enter the interpolation functions (see Equation (12)]. 
One of the four rotational differences is dependent on the other three rotational dif- 
ferences, and hence, there are only three independent rotational degrees of freedom. 
This produces a rank deficiency in the stiffness matrix and leads to a spurious zero- 
energy mode. This spurious mode is associated with a state of zero nodal in-plane 
displacements and equal nodal in-plane rotations. 

However, this mode can be suppressed easily. MacNeal and Harder [16] eliminate 
this mode by adding an energy penalty to the stiffness matrix. The energy penalty is 
defined in terms of the weighted average of the differences between corner rotations 
and the rotation at the center of the element. The rotation at the center of the element 
is obtained by evaluating the skew-symmetric part of the displacement gradient at 
the center of the element. A simpler way to eliminate this spurious mode, but which 
calls for the awareness of any user of the finite element, is to prescribe the value of the 
in-plane rotation at one node in the entire domain being analyzed. This will remove 
the rank deficiency in the membrane part of the stiffness matrix. 

The bending part of the element also has 12 degrees of freedom. Similar to 
the membrane part, 12 independent deformation modes are associated with these 12 
degrees of freedom. Three modes represent rigid-body modes, five modes represent 
constant strain states, the remaining four modes represent higher-order strain states. 
A spurious mode similar to the one in the membrane part does not occur here, though 
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the interpolation function for out-of-plane displacement w [Equation (11)] contains 
out-of-plane (bending) rotational differences. In the membrane part, the rotational 
degrees of freedom, as rotational differences, enter the formulation through the inter- 
polation functions for in-plane displacements and not independently. However, in the 
bending part, the out-of-plane rotations enter the formulation both independently 
[2 nd and 3 rd of Equation (11)] and through the interpolation function for w. Thus, 
the bending part of the stiffness matrix is not rank deficient. 

Assumed-Stress Field 

In general, the selection of stress field approximation is governed by a few 
conditions and guidelines [20] such as 

• the satisfaction of equilibrium equations, 

• interelement traction reciprocity, 

• invariance and completeness of polynomial expansion, and 

• spurious zero-energy modes. 

In earlier hybrid models based on the principle of minimum complementary en- 
ergy [21], the assumed-stress field satisfied equilibrium equations exactly, and conse- 
quently the stress field was expanded in Cartesian coordinates. However, in the hybrid 
formulations such as the present one, the assumed-stress field can be expanded in el- 
ement natural-coordinate basis and the equilibrium equations need not be satisfied 
a priori. In any case, the equilibrium equations are satisfied in a variational sense. 
The invariance of the assumed-stress field could be achieved by having a complete 
polynomial expansion [20]. However, the expansion of assumed-stress field in element 
natural-coordinate basis provides a way to achieve the necessary invariance property 
to the approximations. 

Finally, the chosen stress field should produce no spurious zero-energy modes. Let 
N t be the number of independent stress parameters in the stress field approximations; 
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let N d be the number of deformation modes associated with the element; let N r be 
the number of rigid-body modes associated with the element. In order to suppress 
each independent deformation mode, it is necessary that a corresponding stress term 
exists such that the deformation mode does not produce zero strain energy. As 
such, each independent stress parameter should suppress one non-trivial independent 
deformation mode [20, 22], That is, 

N,> (N d -N r ). (16) 

While this statement sets a lower bound on N t , a definitive statement regarding the 
upper bound on N, cannot be made here. In this context, it is useful to define a ratio 
‘q’ given by 

N, 

a N d -N r ' 

Based on Equation (16), it is clear that a > 1. At this stage, a stress field that is de- 
scribed by a complete polynomial expansion in element natural-coordinate basis and 
that does not satisfy equilibrium equations a priori might appear suitable. However, 
it has been observed that a very high value of a produces a very stiff element [23]. 
Fraejis de Veubeke’s limitation principle [25] states that when a complete indepen- 
dent stress field that does not satisfy equilibrium equations a priori is assumed, the 
hybrid element would reduce to a fully integrated displacement-based element. One 
way to overcome this is to impose additional constraints by enforcing equilibrium 
in the variational statement using Lagrange multipliers [3]. The lower bound on a 
in mixed/hybrid formulations is stated mathematically by the inf - sup condition, 
popularly known as the LBB condition [24]. 

In the present approach the assumed-stress field is not complete, and the equi- 
librium equations are satisfied exactly for regular (rectangular) element shapes. The 
stress field for the membrane part is assumed to be 

Ni = P\ + fat + fan + fan*, 
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Nt] — 02 + Phi + @71 + 
^ = Pi - PiV - Pli ■ 


( 18 ) 


The stress field is expanded in natural coordinates to make the element less sensi- 
tive to mesh distortion. Also note that the stress field is invariant with respect to 
reference natural coordinates. The stress resultants expressed in natural coordinates 
are transformed to Cartesian coordinates by the contravariant tensor transformation 
laws [3]. The transformation rule for the stress field is 


(X 


OLp 

(«0 


J?u.< i) Jfu. i) <)«.-?) 


(19) 


where J * are terms from the Jacobian matrix, given by 


dx n 



and and represent the stress state in Cartesian coordinates and natural 
coordinates, respectively. For a distorted mesh, will be linear in £ and/or 7/ and 
hence the transformation involving jf will have higher-order terms in it. Moreover, 
the constant stress terms will no longer be preserved which causes the element to 
fail the constant strain patch test. In order to preserve the order of stress field 
chosen, including the constant terms, J^’s are evaluated at the origin of the natural- 
coordinate system of the element. Now the transformation rule is 


= J?(0,0) J?(0,0) <({,,). (20) 

This transformation is used in evaluating both the [H] and [ T ] matrices [Equation (5)] 
needed for the linear stiffness matrix [ K } [Equation (8)]. 

The stress field described in Equation (18) will not satisfy equilibrium equations 
a priori for any general, irregular shape. However, this field satisfies equilibrium 
equations a priori for regular (rectangular) element shapes. In any case, as discussed 
earlier, the stress field will satisfy the equilibrium equations in a variational sense. 
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The bending part of the element also has 12 degrees of freedom and hence 12 
independent deformation modes. Three modes are rigid-body modes. Similar to the 
membrane part, 9 independent stress parameters are required to suppress all the other 
deformation modes. Though 9 independent stress parameters are enough to suppress 
all the deformation modes, the following 13-parameter stress field for the bending 
part is less sensitive to mesh distortion, 

= A + A£ + fie 7 ] + 0s 7 ] 2 , 

Mr, = A + A£ + A 7 ? + A£ 2 > 

= A + foot + 0UV + 2^ 13 ^ 2 ’ 

This 13-parameter stress field is a little different from the one used in reference [18]. 
Numerical experimentations indicate that this field is somewhat more accurate than 
the one used in reference [18]. 

The transverse shear stress resultants are obtained by relating the transverse 
shear stress resultants to bending stress resultants, as 

Qt — @4 + 011 + 013 7 ], 

Qr) = 0? + Ao + 012%- (22) 

Note that the transverse shear stress resultant field is coupled to bending stress resul- 
tant field and the equilibrium equations are satisfied a priori for regular (rectangular) 
shapes. The stress fields [Equations (21), (22)] are transformed to Cartesian coordi- 
nate system using the same contravariant tensor transformation laws used for the 
membrane part [Equation (20)]. Also note that a = 1 [Equation (17)] for membrane, 
a = 1.44 for bending, and a = 1.22 for membrane and bending combined. 
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Geometric Stiffness Matrix 


The geometric stiffness or initial stress matrix is derived by using the full nonlinear 
Green-Lagrange strain tensor 

- du k duk ' 

2 dx* 5xj 

Using Reissner-Mindhn kinematic assumptions and retaining the nonlinear terms only 
for the membrane strains, the following relations result. The strains are written as 



{£} = {<«} + *{«} 


( 24 ) 


where {e°} are the membrane strains, {k} are the bending strains, and z is the 
coordinate in the direction through the thickness of the plate or shell. The membrane 
strains are written a s 


in = ui) + u°nl) 


where the linear part is given by 



1 

du°l dx 

{4} = < 

€ v | = < 

dv°j dy 


1 J 

[ du° j dy + dv° j dx 


and the nonlinear part is given by 


{ £ yvx,} ~ 


1 

2< 


/ 8u° \2 

}&* 

v By ) 

0( ® u ° 

Z v ax dy 


+ 

+ 

+ 


/ dv° \2 
v a* / 


+ (Sr ) 2 

_l_ ( Bw c \2 

^ v By f 

(Jv° I 0w° Bw° \ 

Bx 8 m ' 8x 8 m ' 


(dv° \2 
1 By ) 

Bv g Bx 
Bx By 1 Bx By 


The bending strains are computed using the changes in curvature by 


j 

{ ddjdx 

< Ky . 

> = -ddjdy 

, J 

{ dojdy - dejdx , 


The transverse shear strains are given by 


{7“} 


Tj, \ _ / + dw°/9z \ 

-tf, j l -9,+dw°/dy J‘ 


( 25 ) 


( 26 ) 


( 27 ) 


( 28 ) 


( 29 ) 
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The generalized Helhnger-Reissner functional including the nonlinear strains can 


be written as 

Iltf* = j v M t {D]{o} dV + J v {<t) T {C \{ «} dV 

+ L{«o} T {c° kl } dV - / {«} T {t„} dS , (30) 

Jv Js 9 

where {<7 0 } is the prescribed pre-buckling stress state. Upon substitution of the 
displacement and stress approximations, Equation (30) reduces to 

+ larmw - wqn + oi) 

where the element geometric stiffness matrix [ K „.] is given by 

IK'} = Jim T {G} T {SmiN} dA (32) 

J A 

and 

'■a, 0 0 0 0 0 ' 

d y 0 0 0 0 0 

fr] _ o d x o o o o 

. l J 0 dy 0 0 0 0 

0 0 d x 0 0 0 

.0 0 dy 0 0 0 . 

where d x = d/dx and d y — d/dy. The matrix [«S], which corresponds to the pre- 
buckling stress state, consists of membrane stress resultants which are either pre- 
scribed or evaluated from a linear static stress analysis. It is given by 

' S 0 0 ' 

[5] = 0 S 0 (33) 

0 0 s _ 

where 

N° N° 1 

O _ 1, x 1J xy 

N° N° 

*v v\ 

Consistent Mass Matrix 


In linear free vibration analysis, it is assumed that the system is undamped and there 
are no external forces. The Hellinger-Reissner functional [Equation (1)] with the 
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kinetic energy contribution can be written as 


n " B = C JyWMW dV + (W'WW' 11 ' 

- \ j v {u} T [p]{i) iv) it, (34) 

where dots (.) over u indicate time derivatives. Upon substitution of stress field and 
displacement field approximations, Equation (34) reduces to 

n « B = + w r m{?} - <« , ( 35 ) 

where the consistent mass matrix [M c ], representing the inertia of the system, is given 
by 

[JWJ = ! [N\ T {P\{N] iA (36) 

J A 


and 


W5X6 — 


m 0 0 0 0 m.i 

0 m 0 0 —m x 0 

0 0 m 0 0 0 

0 — mi 0 m 2 0 

m-i 0 0 0 m 2 


(37) 


where 


m x 


m 2 


• - 4 

- 4 


p dz , 


pz dz , 


= J' t pz’dz 


(38) 


and p is the mass density. 


Numerical Results 


The element described herein was developed and implemented within the framework 
of NASA Langley Research Center’s computational structural mechanics testbed 
COMET [26] using the generic element processor template [27]. The element will 
be referred to herein as A4S1. 
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The proper selection of displacement and stress fields has automatically elimi- 
nated spurious zero-energy modes and locking in the present element and has made 
the element invariant with respect to reference coordinates. However, the element 
still needs to be validated with regard to its convergence and its sensitivity to mesh 
distortion. MacNeal and Harder [28] have collected and proposed a standard set of 
problems to test and validate a new element. The element’s performance is observed 
by analyzing these problems and a few other classical problems. The results are com- 
pared with known theoretical solutions and other numerical solutions. Other 4-node 
quadrilateral shell elements used for comparison are briefly described in Appendix A. 
Consistent units have been used throughout the analyses. 

Patch Test 

This test was originally proposed by Irons [13] to establish the convergence of 
an element. The patch test as shown in Figure 1 involves a specific mesh of distorted 
elements with specified applied displacements that result in a constant state of either a 
membrane or bending stress field. It is now widely recognized that an element should 
pass the patch test, failing which, its convergence will not be assured. A4S1 passes 
both the membrane and bending patch tests. The recovered strains and stresses are 
exact. 

MacNeal-Harder Straight Cantilever Beam 

This test case was proposed by MacNeal and Harder [28]. The capability of the 
element to handle constant and linearly varying strains and curvatures is tested by 
applying appropriate unit loads at the free end. The beam is modeled by a 6 x 1 
mesh. It is discretized using rectangular-, trapezoidal- and parallelogram-shaped 
elements. The element’s sensitivity to mesh distortion is studied and compared with 
other elements. The finite element discretizations and the exact solutions [28] are 
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shown in Figure 2. The theoretical result for twist (load D) is 0.03406 as reported 
in Timoshenko and Goodier [29]. However, reference [28] reports the solution to 
be 0.03028. Analysis with successively refined meshes converged to 0.03385 which 
is much closer to Timoshenko and Goodier solution [29], and hence this solution is 
used here for normalization. The normalized results given in Table I indicate that 
the elements predict good results for rectangular-shaped elements and the results 
deteriorate for the two distorted mesh cases. However, A4S1 performs very well and 
can handle distorted mesh configurations effectively. 

• Remark 1: Attention is drawn to the case of in-plane shearing load (load B), 
particularly in distorted mesh configurations. 4.ANS, 4.HYB, and QUAD4 
perform badly. The advantage of drilling degrees of freedom is evident from 
the superior performance of the present element, the performance of Q4S, and 
a relatively modest improvement shown by 4_STG. 

• Remark 2: A4S1 has shown improvement over AQR8 and AQD4 for the twist 
(load D), and this is attributed to the modified bending stress field. 

• Remark 3: The rotational degrees of freedom by virtue of their presence in 
the element displacement shape functions, enter into the calculations of work- 
equivalent nodal loads. Hence, the in-plane and out-of-plane loads at the beam 
tip correspond to both statically equivalent loads and the work-equivalent loads 
including moments corresponding to the rotational degrees of freedom. 


Twisted Cantilever Beam 

This test case involves a straight cantilever beam twisted 90° over its length [28] . 
The twisted beam shown in Figure 3 is modeled by 12 elements along the length and 
2 elements across the width. Each element has a 7.5° warp and the effect of this warp 
on the element’s performance is studied. The results tabulated in Table II indicate 
that A4S1 performs very well, but does not show improvement over AQR8 and AQD4. 
While 4.STG and QUAD4 perform well, 4_ANS and 4.HYB do not perform very well. 
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Curved Cantilever Beam 


The curved cantilever beam formed by a 90° circular arc is shown in Figure 4. 
The loading consists of in-plane shear and out-of-plane shear at the tip [28]. The 
results tabulated in Table III indicate that the new element performs very well. Again 
A4S1 has shown a modest improvement over AQR8 and AQD4 for the out-of-plane 
shear (load C) and this is again attributed to the improved bending stress field. 

Scordelis-Lo Roof 

The Scordelis-Lo roof [30] is a singly curved shell structure (see Figure 5). 
The shell has a gravity load or a uniform dead load (UDL), and the parameter of 
interest is the vertical displacement at the mid-point of the free edge (point A in 
Figure 5). Though the theoretical value used in reference [30] is 0.3086, the value 
used by MacNeal and Harder [28] is 0.3024. This latter value is used for normalization 
here. Due to the symmetry of the structure and loading, only one quadrant of the shell 
is modeled. The convergence behavior is studied for meshes from 2 x 2 to 10 x 10 and 
is tabulated in Table IV. The A4S1 element shows a rapid and monotonic convergence. 

Morley’s Spherical Shell 

The spherical shell is doubly curved (see Figure 6) and the equator of the shell is 
chosen to be a free edge. Hence, the problem reduces to a hemisphere with four point 
loads alternating in sign at 90° intervals along the equator. An 18° hole has been 
introduced at the top of the hemisphere to avoid having to model the pole [28]. The 
theoretical solution for displacement at the point of load and in the direction of the 
load is reported as 0.094 in reference [28] and it is used here for normalization. Only 
a quadrant of the shell is modeled. Convergence is studied for meshes from 2 x 2 to 
12 x 12 and the results are tabulated in Table V. While 4.ANS and QUAD4 seem 
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to show a very rapid convergence, A4S1, AQR8, AQD4, and 4.STG exhibit a slow 
convergence. Note that even for a 12 x 12 mesh, A4S1 shows an 8.4% error. 

• Remark 4: The present problem involves significant contribution from both 
membrane and bending strains to the radial displacement at the points of 
loading. The reason attributed for the slow convergence of the elements with 
drilling rotations [16] is that the stiffness of the shell is increased by membrane- 
bending coupling arising out of the faceted geometric approximation of the 
shell. This coupling is actually the coupling between drilling rotations and 
bending rotations due to the changes in slopes at element interfaces (because 
of faceted representation). 

• Remark 5: Furthermore, the membrane-bending coupling may be amplified by 
the fact that the drilling rotations are not true rotations, but are rotational 
connectors. An independent drilling rotation entering the formulation either 
through the variational statement or through the kinematics may increase the 
flexibility of the discretized shell. 


Isotropic Square Plate 

The classical plate bending problem of a clamped, isotropic square plate sub- 
jected to a concentrated load at its center is analyzed. The convergence behavior, 
the effect of mesh distortion, and the effect of varying the thickness of the plate are 
studied. The geometry of the plate and its properties are shown in Figure 7. The 
maximum deflection at the center of the plate, as calculated in reference [31] using 
the Ritz approximation, is given by 

= 0.0056^ (39) 

where D is the flexural rigidity of the plate, P is the applied concentrated load for 
the full plate, and a is the side length of the plate. The central deflections obtained 
from the finite element analyses are normalized with respect to the value calculated 
from Equation 39. 

The convergence behavior shown in Figure 7 indicates that A4S1 converges 
rapidly and monotonically. Next, the analysis is carried out using a distorted 7x7 
mesh. The center of the quarter plate is fixed and the nodes immediately adjacent to 
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the central node are rotated from 0 to 45 degrees. The boundary nodes are also fixed. 
From Figure 8, the mesh distortion is shown to have the least effect on A4S1. Again 
A4S1 offers a modest improvement over AQR8 and AQD4. As mentioned earlier, the 
consistent interpolation order between the transverse displacement and out-of-plane 
rotations eliminates shear locking in A4S1. To confirm this with numerical valida- 
tion, the effect of decreasing the thickness of the plate on the element’s performance 
is studied. The length-to-thickness ratio (defined as a/h ) is varied from 10 to 10,000. 
None of the elements lock at high length-to-thickness ratio. Representative results 
for two elements are shown in Figure 9. 

Pear-Shaped Cylinder 

This is an extremely interesting problem because of its continually varying ge- 
ometry. The cross-section of the cylinder is shown in Figure 10, which supposedly is a 
representation of an early space shuttle fuselage configuration. The shell is isotropic 
with a uniform thickness of 0.01 inches. Only one quarter of the cylinder is modeled 
due to the symmetry of the structure. Symmetry boundary conditions are imposed 
on three edges and a uniform end shortening is applied to the simply supported edge. 
The analysis is carried out using three different meshes: 4 x 23 (92 nodes, 66 ele- 
ments), 5 x 34 (170 nodes, 132 elements), and 7 x 51 (357 nodes, 300 elements), 
where n x m refers to n nodes along the length of the cylinder and m nodes along 
the half-cylinder circumference. However, the solution reported here is for the 7x51 
mesh, where the solutions converge. 

One of the earliest available numerical analyses on this problem was performed 
by Hartung and Ball [32]. The variation of radial displacements and the axial stress 
resultants with 6 measured counter-clockwise from the top, are shown in Figures 11 
and 12, respectively. Results of 4_ANS and A4S1 are very close. Further comments 
on the accuracy of the solution cannot be made here due to the lack of Unear static 
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analysis data on this problem. However, this problem has been an attractive target 
for nonlinear analyses and elastic shell collapse studies, due to its complex nonlinear 
collapse behavior. Herein a linear elastic buckling analysis is carried out on the 
cylinder. The buckling loads for two modes are presented in Table VI and the buckling 
mode shapes are shown in Figure 13. Since there are no analytical solutions available 
for comparison, the results are normalized with respect to converged solutions of a 
9-node assumed natural strain element (9_ANS) [26, 33]. The present element and 
4.STG which contain the drilling degree of freedom, predict lower buckling loads 
compared to 4.ANS. 

Rectangular Plate 

This problem is analyzed to validate the element’s capability to handle linear 
buckling and free vibration problems. This is a classical problem for which theo- 
retical solutions exist. The rectangular plate [see Figure 14] is isotropic and simply 
supported on all sides. The plate is subjected to a uniform uniaxial compression and 
a biaxial compression [see Figure 14(a) and 14(b), respectively]. Only a quarter of 
the plate is discretized, with a 5 x 8 mesh [shaded region in Figure 14(a) and 14(b)]. 
The linear buckling analysis is performed with the pre-buckling stress state consist- 
ing of stresses computed from a linear static analysis. The results are tabulated in 
Table VII and VIII, where m and n refer to the number of half waves in X and Y 
directions, respectively. The results are normalized with respect to the theoretical 
solutions reported in reference [34]. A4S1 predicts quite accurate buckling loads and 
again, A4S1 and 4.STG predict lower buckling loads than 4_ANS. The plate is also 
subjected to an uniform shear and the full plate is discretized with a 5 x 8 mesh [see 
Figure 14(c)]. The results are shown in Table IX. A4S1 predicts a lower and more 
accurate buckling load than 4.ANS. 
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A free vibration analysis is carried out on the plate in its unstressed state. The 
results are tabulated in Table X, where m and n refer to the number of half waves 
in X and Y directions, respectively. The results are normalized with respect to the 
theoretical solutions reported in reference [34]. A4S1 gives quite accurate natural 
frequencies, but 4_HYB seems to have performed marginally better. 

Axially Compressed Cylinder 

A cylindrical shell with simply supported edges is subjected to an uniform axial 
compression (see Figure 15). The prestress state is computed from a linear stress 
analysis of the axially compressed cylinder. Note that this prestress state is not 
exactly constant due to the simply supported boundaries, which prevents the shell 
from uniformly expanding in circumferential direction. Another interesting feature 
of this problem is the closely packed eigenvalues and the corresponding considerably 
different mode shapes. Modeling the entire cylinder is computationally intensive and 
hence only a 15-degree sector and one tenth of the length of the cylinder, sufficient to 
capture the lowest mode, is modeled. Symmetric boundary conditions are employed 
on all edges except the loaded edge which is simply supported. A convergence study 
is done for meshes N x N, where N (number of nodes per side) ranges from 3 to 9 
(see Figure 15). A4S1 and 4.STG show a better convergence than 4_ANS. The first 
two mode shapes of this discretization are shown in Figure 16. 

Cylindrical Shell 

An octant of a cylindrical shell (see Figure 17) is analyzed for the free vibrational 
characteristics of the cylindrical shell. To capture the correct modes, symmetric and 
anti-symmetric boundary conditions are defined for the edge BC. Edge AB is simply 
supported and Edges AD and DC are in the symmetric planes. The results are 
compared with the theoretical solutions reported in reference [35] and are tabulated 
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in Table XI. The mesh refers to number of nodes in the circumferential direction by 
number of nodes in the axial direction. A4S1 and 4_HYB converge to the theoretical 
solution. The vibration mode shapes are shown in Figure 17. 

Conclusions 

Numerical results indicate that the element developed in this paper is devoid of the 
usual deficiencies of 4-node shell elements. Though the element has one spurious 
zero-energy mode, this mode can be suppressed easily by prescribing the value of 
in-plane normal rotation at one node in the entire finite element model. It has been 
shown that the element is nearly insensitive to mesh distortion, does not lock, and 
has desirable convergence and invariance properties. The element has also performed 
very well for buckling and free vibration analyses. 
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Appendix A 


EX47(ES1): A C°, isoparametric, assumed natural-coordinate strain element, devel- 
oped by Park and Stanley [33]. The element is not invariant with respect to 
reference coordinates and does not pass the patch test. This element will be 
referred to as 4_ANS. 

QUAD4: A 4-node isoparametric shell element with selective reduced integration [28]. 
This element is available in MSC/NASTRAN. Q4S, an improved version of 
QUAD4, contains the drilling degree of freedom [16]. 

E410(ES5): A C 1 , incompatible, displacement based element with drilling degrees of 
freedom. The element kernels were extracted from the STAGS computer code 
(Almroth et al.[36]). The element uses cubic interpolation for the displacement 
field. This element is not invariant and does not pass the patch test. This 
element will be referred to as 4_STG. 

EX43(ES4): A ( 7 °, isoparametric^ assumed-stress hybrid shell element with no drilling 
degree of freedom. This element is invariant with respect to reference coordi- 
nates and passes the patch test. This element was developed by Aminpour[37]. 
This element will be referred to as 4_HYB. 

AQR8(ES8): An assumed- stress hybrid shell element with drilling degrees of freedom, 
developed by Aminpour [18]. The formulation is based on Hellinger-Reissner 
variational principle. The 4-node element is obtained from an “internal” 8- 
node element. 

AQD4(ES8): An assumed-stress hybrid shell element with drilling degrees of free- 
dom, developed by Aminpour [19]. The formulation is based on modified 
complementary energy principle and the stress field is expanded in Cartesian 
coordinates. 
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Table I. MacNeal-Harder straight cantilever beams. 


A : Extension, B : In-plane Shear, C : Out-of-plane Shear, D 

: Twist 

Load 

4.ANS 

QUAD4 

4.STG 

4.HYB 

AQR8 

AQD4 

A4S1 

Rectangular-Shaped Elements 

A 

0.995 

mzm 

0.994 

0.996 

m 

0.998 

0.998 

B 

0.904 

0.904* 

0.915 

0.993 

Hi 

0.993 

0.993 

C 

0.980 

0.986 

0.986 

0.981 

0.981 

0.981 

0.981 

D 

0.856 

0.941** 

0.680 

1.023 

1.011 

1.011 

1.009 

Trapezoidal-Shaped Elements 

A 

WHOM 

0.996 

0.991 

0.999 

0.998 

0.998 

0.998 

B 

IBS 

0.071* 

0.813 

0.052 

0.986 

0.986 

0.986 

C 

BEH 

0.968 

0.925 

0.075 

0.965 

0.965 

0.969 

D 

0.843 

0.951** 

0.683 

1.034 

1.029 

1.009 

1.007 

Paralle 

ogram-S 

raped Elements 

A 

0.966 

EEZB 

0.989 

0.999 

0.998 


sm 

B 

0.324 

0.080* 

0.794 

0.632 

0.977 


Hi 

C 

0.939 

0.977 

0.991 

0.634 

0.980 


0.980 

D 

0.798 

0.945** 

0.677 

1.166 

1.159 

1.010 

1.007 


* Q4S results for these cases are 0.993, 0.988, and 0.986, respectively. 

** These results were normalized with respect to 0.03208 in reference [28]. How- 
ever, all the other results for twist are normalized with respect to 0.03406 as 
reported in Timoshenko and Goodier [29]. 


Table II. MacNeal-Harder twisted cantilever beam. 


Load 

4.ANS 

QUAD4 

4.STG 

4.HYB 

AQR8 

AQD4 

A4S1 

In-plane shear 
Out-of-plane shear 

1.357 

1.293 

■1 

1.054 

1.173 


m 

m 

0.991 

1.093 


Table III. MacNeal-Harder curved cantilever beam. 


Load 

4.ANS 

QUAD4 

4.STG 

4.HYB 

AQR8 

AQD4 

A4S1 

In-plane shear 
Out-of-plane shear 

m 

m 

B 

m 

H 

B 

0.997 

0.970 




























































Table IV. Scordelis-Lo roof. 


Mesh 

4_ANS 

QUAD4 

4.STG 

4_HYB 

AQR8 

AQD4 

A4S1 

2x2 

1.387 

1.376 

1.384 

1.459 

1.218 

1.218 


4x4 


1.050 



1.021 

1.021 


6x6 

■E9 

1.018 

1.015 

1.028 

1.006 

1.006 

1.007 

8x8 

1.005 

1.008 

1.005 

1.017 

1.003 

1.003 

1.003 

10 x 10 

1.003 

1.004 

1.001 

1.011 

1.001 

1.001 

1.001 


Table V. Morley’s Spherical Shell. 


Mesh 

4.ANS 

QUAD4 

4.STG 

4_HYB 

AQR8 

AQD4 

A4S1 

2x2 

0.968 

0.972 

0.338 


0.382 

0.381 

| Ml 

4x4 

1.018 

1.024 

0.519 


0.227 

0.226 

ms 

6x6 

1.001 

1.013 

0.841 

1.060 

0.432 

0.432 

0.433 

8x8 

0.995 

1.005 

0.949 

1.040 

0.681 

0.680 

0.682 

10 x 10 

0.993 

1.001 

0.978 

1.027 

0.835 

0.835 

0.838 

12 x 12 

0.992 

0.998 

0.988 

1.020 

0.914 

0.914 

0.916 


Table VI. Linear buckling analysis of the pear-shaped cylinder. 




^cr / ^cr 

Modes 

9-ANS 





— raa 

4.393 

USM 

nun 

MEjm 

|R|H| 

2 

6.483 

mi 

Qy| 


BUS 


Table VII. Linear buckling analysis of the rectangular plate under uniaxial 

compression. 




N* 

V 

", / K 

Modes 

71 

Analytical 

4_ANS 

4_STG 

4.HYB 

A4S1 

1 

1 

523.04 


0.991 

■■•Mil 

0.998 

2 

3 

874.89 


0.988 

mrn&m 

1.001 














































Table VIII. Linear buckling analysis of the rectangular plate under biaxial 

compression ( N x = N y = N). 


Modes 

m 

n 

N* 

Analytical 

N / N* 

4.ANS 

4.STG 

4_HYB 

A4S1 

1 

i 

i 


| 

■Egg 

mm 

0.998 

2 

i 

3 

| 



WBm 

1.003 

3 

3 

1 

1152.36 



HU 

1.020 


Table IX. Linear buckling analysis of the rectangular plate under uniform shear. 


Modes 

N* 

xy 

Analytical 


4_ANS 4_STG 4JIYB A4S1 

1 

3629.69 

1.284 0.930 1.213 0.989 


Table X. Linear free vibration analysis of the rectangular plate. 


Modes 

m 

n 

Frequency 

r 

Analytical 

fir 

4.ANS 

4.HYB 

A4S1 

1 

i 

i 

112.40 

| 

■ 

IMI 

2 

i 

3 

436.10 

1 

HR® 

1.037 

3 

3 

1 

687.80 

1.111 


■flKM 


Table XI. Linear free vibration analysis of the cylindrical shell. 


Number of circumferential waves = 6 


Mesh 

Symmetric Mode 1 

Anti-symmetric Mode 1 


HnSgH 


u> 2 /lJT 

4.ANS 

4.HYB 

A4S1 

4.ANS 

4.HYB 

A4S1 

25 x 9 
31 x 11 
35 x 13 

1.022 

1.008 

1.002 

H 

mu 

m 

HH 

■ 

■ 

1 - 
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Figure 2. MacNeal- Harder straight cantilever beams. 
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12.0 

Width 

1.1 

Depth 

0.32 

Elastic Modulus 

29.0 x 10 6 

Poisson’s Ratio ( i /) 

0.22 



Figure 3. MacNeal-Harder twisted cantilever beam. 








Inner Radius 


Outer Radius 


Depth 


Elastic Modulus 

10 7 

Poisson’s Ratio (u) 

0.25 



Figure 4. MacNeal-Harder curved cantilever beam. 








L = 50, R = 25, t = 0.25 
E - 4.32 x 10 8 , v = 0.0 
UDL = 90, u = w = 0 on curved edges 



Figure 5. S cor delis- Lo roof. 


Radius = 10., Thickness = 0.04 
E = 6.825 x 10 7 , v = 0.3 



Figure 6. Morley’s spherical shell. 
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Figure 7. Convergence study for the isotropic square plate; a = 4, h = 0.001, 
P = 100, E = 10 7 , ^ = 0.3. 
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Figure 9. The effect of thickness change (for a 7 x 7 mesh). 



Figure 10. The pear-shaped cylinder. 
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Mode 1 

N y / N y (converged) = 1.013 



Mode 2 

N y / N y (converged) = 1.016 


Figure 13. Buckling mode shapes for the pear-shaped cylinder. 




Figure 14. Linear buckling analysis of the rectangular plate under (a) uniaxial 
compression, (b) biaxial compression, and (c) uniform shear; length along the 
X-direction = 15, length along the Y-direction = 20, E = 30 x 10 6 , v = 0.3. 






Discretized Model (9 x 9) 




Mode 1 
>»cr = 0.99892 


Mode 2 
A-cr = 1 .02980 


Figure 16. Mode shapes corresponding to the lowest two critical loads 
for the axially compressed cylinder. 
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Symmetric Mode 1 
to 2 = 0.636 x 10 5 


Anti-symmetric Mode 1 
to 2 = 0.584 x 10 5 


Figure 17. Vibration mode shapes for the cylindrical shell; E = 10 7 , v - 0.3, 
p = 0.1, R = 10, t = 0.1, h = 5. 



